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Due to lack of scientific understanding, some mechanisms may be missing in 
mathematical modeling of complex phenomena in science and engineering. These 
mathematical models thus contain some uncertainties such as uncertain parame- 
ters. One method to estimate these parameters is based on pathwise observations, 
i.e., quantifying model uncertainty in the space of sample paths for system evo- 
lution. Another method is devised here to estimate uncertain parameters, or 
unknown system functions, based on experimental observations of probability 
distributions for system evolution. This is called the quantification of model un- 
certainties in the space of probability measures. A few examples are presented to 
demonstrate this method, analytically or numerically. 
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1 . Introduction 

In this chapter we discuss some issues about quantification of model uncertainties 
in complex dynamical systems. 
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Mathematical models for scientific and engineering systems often involve with 
some uncertainties. We may roughly classify such uncertainties into two kinds. The 
first kind of uncertainties may be called model uncertainty. They arc due to physical 
processes that are not well understood or not well-observed, and thus are not or 
not well represented in the mathematical models. 

The second kind of uncertainties may be called simulation uncertainty. This 
arises in numerical simulations of multiscale systems that display a wide range of 
spatial and temporal scales, with no clear scale separation. Due to the limitations of 
computer power, not all scales of variability can be explicitly simulated or resolved. 

These uncertainties are sometimes also called unresolved scales, as they are not 
represented (i.e., not resolved) in modeling or simulation. Although these unresolved 
scales may be very small or very fast, their long time impact on the resolved sim- 
ulation may be delicate (i.e., may be negligible or may have significant effects, ^^21121 
or in other words, uncertain). Thus, to take the effects of unresolved scales on the 
resolved scales into account, representations or parameterizations of these effects 
are desirable. 

Model uncertainties have been considered in, for example,. ^i^MIillilHlll Research 
works relevant for parameterizing unresolved scales include,! ^^*^^ * ^^ * ^^ * '^^ * '^^ * ^" * ^^ ^^ 
among others. Stochastically representing unresolved scales in fluid dynamics has 
considered as well.'^SBnill! 

In this chapter, we only consider model uncertainties. Specifically, we consider 
dynamical systems containing uncertain parameters or unknown system functions, 
and examine how to estimate these parameters, using observed probability distri- 
butions of the system evolution. 

After briefly comment on estimating uncertain parameters based on observed 
sample paths for the system evolution in ^ }2~1 we then, in ! }3~1 propose a method of 
estimating uncertain parameters based on observed probability distributions (i.e., 
probability measures) and present a few examples to demonstrate this method, 
analytically or numerically. 

2 . Quantifying uncertainty in the space of paths 

Since random fluctuations are common in the real world, mathematical models 
for complex systems are often subject to uncertainties, such as fluctuating forces, 
uncertain parameters, or random boundary conditionsi ^^ l ^^ l ^^ l ^'^ l ^^ l ^^ l Stochastic 
differential equations (SDEs) such as 

dX ^b{X)dt + a{X)dBt, (2.1) 

are appropriate models for many of these systems .'^'^^'^ Here Bt is a Brownian 
motion or Wiener process, the drift b{X) and diffusion a{X) contain uncertain 
parameters (or (6(-) and (t(-) are unknown), to be estimated based on observations. 



For example, the Langevin type models are stochastic differential equations de- 
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scribing various phenomena in physics, biology, and other fields. SDEs are used to 
model various price processes, exchange rates, and interest rates, among others, in 
finance. Noises in these SDEs may be modeled as a generalized time derivative of 
some distinguished stochastic processes, such as Brownian motion (BM) or other 
processes. 

We are interested in estimating parameters contained in the stochastic differen- 
tial equation p .ip , so that we obtain computational models useful for investigating 
complex dynamics under uncertainty. 

Theoretical results on parameter estimations for SDEs driven by Brownian mo- 
tion are relatively well developed j "^ * ^ * ^'^ * ^'^ ' ^'^ '^"^ "^-^ and various numerical simulations 
for these parameter estimationJ^'^EIIlTl are also implemented. Sed^ for a more re- 
cent review about estimating and computing uncertain parameters, when dynamical 
systems are submit to colored or no n- Gaussian noises. 

These research works on estimating uncertain parameters in dynamical systems 
are based on observations of sample paths. In the next section, we devise a method 
to estimate uncertain parameters based on observations of probability distributions 
of the system evolution. 

3 . Quantifying uncertainty in the space of probability measures 

Consider a dynamical system with model uncertainty, modeled by a scalar SDE 

dX = b{X)dt + a{X)dBt, X{0)=xo, (3.1) 

where the drift b{X) and diffusion a{X) contain uncertain parameters, to be esti- 
mated based on observations of probability distributions (i.e., probability measures) 
of the system paths Xf . 

To this end, we need to introduce the Hellinger distanc^^^ between two prob- 
ability measures. It is used to quantify the similarity between two probability 
distributions. This is a metric in the space of probability measures. 

For our purpose here, we define the Hellinger distance H{f,g) between two 
probability density functions p{x) and q{x) as follows 

H'iP, 9) " ^ / iVm- ?dx. (3 .2) 

The Hellinger distance H{p,q) satisfies the property: < H{p,q) < 1. 

We estimate uncertain parameters by minimizing the Hellinger distance between 
the true probability density p for the solution process X(t) and its observed proba- 
bility density q. In reality, the probability density p has to be numerically formulated 
or discretized. But in order to demonstrate the method, we consider two examples 
for which the true probability density p can be analytically formulated. In the first 
example, we minimize the Hellinger distance between the true stationary proba- 
bility density for the solution process X{t) and its observed stationary probability 
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density, while in the second example, we do this for time-dependent probability 
densities. 

3.1. Observation of stationary probability distributions 

Under appropriate conditions on b and a {see^^ p. 170), such as, 6 < and a ^ 
as well as some smoothness requirements, there exists a stationary probability 
density p{x) for the SDE p .ip . as a solution of the steady Fokker-Planck equation 
Pa:x + {hs.m{x)p)x = 0, 

C r, ^dy 

where the positive normalization constant C is chosen so that p > and J^p(x)dx = 
1, i.e., 

C^ll , dx. 

Note that x* here may be an arbitrary point so that the integral j^, ^^dy exists 
(say, take x* = if that is possible). 

Example 3 .1. (i) A special case: Langevin equation 



dX = ~bXdt + dBt, 

with parameter & > 0. Given an "observation" of the stationary probability density 
q[x) = -^e~^ . Find a 6 so that the Hellinger distance F{b) = | J^[\/pix) — 

^ qixy^dx is minimized, 
(ii) A more general case: 



dX = b{X)dt + dBt, 

with function b{x) < 0. Given an "observation" of the stationary probability density 
q{x) = ^ ^^^i (the Cauchy distribution). Find a function b{x) < so that the 
Hellinger distance F{b{x)) = i J-^[\/ p{x) ~ \/ q{x)]'^dx is minimized. 

Solution: 

(i) The true stationary probability density for the solution process Xt is 

Insert p, q into the Hellinger distance F{b), which is now an algebraic function of 
parameter 6 > 0. Thus wc use deterministic calculus to find a minimizer b (possibly 
by hand, or Matlab if needed). Note: /^e"^ dz = -y/7r. 
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To minimize the Hellinger distance F(b), we calculate its derivative 

^ Jr 2Vnb VTT VT^ JR 4 2 

= -i=6^(i + fe)-i__l=6-|(i + 6)-^ ^0. 

Therefore, 

64(1 + 5) 2 = . 

Thus we obtain the parameter 6=1. 

(ii) The true stationary probability density for the solution process Xt is 

g2 £ b(y)dy 



p{x) 



r°° ^2f;^b{y)dydx- 
J — oo 



Insert p,q into the Hellinger distance F{b{x)), which is now a functional of 
b{x) and thus we use calculus of variations (on F{b{x))) to find a minimizcr b{x). 
We then derive the Euler-Lagrangc equation to be satisfied by 6(x), together with 
appropriate boundary conditions (needed for p{x) > and /jgp(a;) dx = 1). 

To this end, we calculate, for an arbitrary "variations" h{x) 



F{b{x) + eh{x)) 



(,'^IoiKv)+eh{y))dy 2e^o i^iy)+^Hv))dx 

dx. 



J^^^2j„-ib(y)+eHy))dy^^ ^^-—^^^^^ /o^ (;>(y)+eMy))<i,^^ 7r(l 4 

The Euler-Lagrange equation for b{x) comes from: ^Ie^o F{b{x) + eh{x)) = 
for arbitrary "variations" h{x). In fact, the Euler-Lagrange equation for b{x) is 

efS ^fo _ V-^« . / (e^Io K-)''- / h{y)dy] dz = 

\ ^7r(l + a;2) j Jr\ Jz ) 

After changing the order of integration (first on y and z then on y and x), we 
have 




rr~^^jp^(w)d^^\ 

fo b{^o}du. I g/o^ bMdu, _ V \dx - I /o K'^)d^dz \ h[y)dy = 

^7r(l + a;2) 



holds for all /i(?/). 
Therefore, 
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Since e^/o K™)^^^^ > q, we further obtain 

g/J= b{w)dw b(w)dw _ \1 \dx = Q 

for y e _R. Then taking the derivative of the above equation with respect to we 
arrive at 



Thus, after taking the 'square' and ' In' on both sides of the above equation, we 

get 



2 / b{w)dw^\n{ e'^J'o''^'"'>''"'dx)-\n{TT{l+y'^j). 







FinaUy taking the derivative with respect to y, we have 

7r(l + y^) 

Also note that we only need b{w)dw < for all y £ B} for the stationary 
probability density to make sense. 



3 .2. Observation of time- dependent probability distributions 

Consider a scalar SDE 

dX = b{X)dt + a{X)dBt, X{0) = xo. (3.3) 

The Fokker-Planck equatiorPSI^IIlll for the probability density p{x, t) = p{x, t; xq,0) 
for the solution X{t, xq) is 

Pt ^ ^{a'^{x)p{x,t))xx - {b{x)p{x,t))^, p{x,0) = S{xo). (3.4) 

With an observation oi p{x,t), we can estimate parameters, or 6(-), or (t(-), by 
examining the inverse problem of the Fokker-Planck equation p .4^ . For more 
information about inverse problems of partial differential equations, see-f^Sl 
Let us look at a specific example. 

Example 3 .2. Consider a scalar SDE 

dX = -b sm{X)dt + V2dBu X(0) = 0. (3 .5) 
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(i) Assume that an observation obtained for p to be 

Find the parameter b by minimizing the Hellinger distance H(jj, qi). 
(ii) Assume that another observation obtained for p to be 

Find the parameter b by minimizing the Helhnger distance H(jj, q2). 
Solution: 

The Fokker-Planck equation for (|3 .5p is 

Pt =Pxx + {bsm{x)p)j;, p{x,0) = S{0). (3 .8) 

In this case, we define the Hehinger distance: 



Hf{b)= max / (Vg^OM) - v^KM7^)^rfa; 

where i = 1,2, and T is the time period when qi{x,t) are observed. We numerical 
find 6 by minimizing Hi{b). 

The observation qi{x,t) is plotted in Figure lUTTl 
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Fig. 0.1. Observation qi(x,t): t = 5 (top) and t = 30 (bottom). 



By the definition Hl{b) — maXigjo^T] /^(V (li{x, t) — yjp{x, t, b))'^dx, we have 
the plot of Hi (b) in Figure 10.21 And whatever T is , Hi (6) is always minimized 
when 6 = 0. This gives us the parameter value 6 = 0. 

The observation q2{x,t) is plotted in Figure 10751 

By the definition -ff|(6) = ^'<^y^tG[o,T] J^ao^y/ 'l2{x , t) — y^p{x, t, b))'^dx, we have 
the plot of H2{b) in Figure [074l So we see that ii T ~ 5, H2{b) is minimized when 
b = 0.7 and if T = 30, H2{b) is minimized when b ~ 0.6. 
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Mathematical models for scientific and engineering systems often involve with 
some uncertainties. We may roughly classify such uncertainties into two kinds. The 
first kind of uncertainties may be called model uncertainty. They arc due to physical 
processes that are not well understood or not well-observed, and thus are not or 
not well represented in the mathematical models. 

The second kind of uncertainties may be called simulation uncertainty. This 
arises in numerical simulations of multiscale systems that display a wide range of 
spatial and temporal scales, with no clear scale separation. Due to the limitations of 
computer power, not all scales of variability can be explicitly simulated or resolved. 

These uncertainties are sometimes also called unresolved scales, as they are not 
represented (i.e., not resolved) in modeling or simulation. Although these unre- 
solved scales may be very small or very fast, their long time impact on the resolved 
simulation may be delicate (i.e., may be negligible or may have significant effects, 
or in other words, uncertain). Thus, to take the effects of unresolved scales on the 
resolved scales into account, representations or parameterizations of these effects 
are desirable. 

Model uncertainties have been considered in, for example,. ^^^I^IIiniilHlll Research 
works relevant for parameterizing unresolved scales include, ^^ * ^^ * ^^ * ^^ * '^" * '^^ * "^^ * ^^ * ^'' ! 
among others. Stochastically representing unresolved scales in fluid dynamics has 
considered as weU.'^SIMII] 

In this chapter, we only consider model uncertainties. Specifically, we consider 
dynamical systems containing uncertain parameters and examine how to estimate 
these parameters, using observed probability distributions of the system evolution. 

After briefly comment on estimating uncertain parameters based on observed 
sample paths for the system evolution in i }2~l we then, in i }3~l propose a method of 
estimating uncertain parameters based on observed probability distributions (i.e., 
probability measures) and present a few examples to demonstrate this method, 
analytically or numerically. 

2 . Quantifying uncertainty in the space of paths 

Since random fluctuations are common in the real world, mathematical models 
for complex systems are often subject to uncertainties, such as fluctuating forces, 
uncertain parameters, or random boundary conditionsi ^^ ' ^^ l ^° l ^^ l ^^ l ^^ l Stochastic 
differential equations (SDEs) such as 

dX ^b{X)dt + a{X)dBt, (2.1) 

arc appropriate models for many of these systems .'^'^^'^ Here Bt is a Brownian 
motion or Wiener process, the drift b{X) and diffusion cr{X) contain uncertain 
parameters, to be estimated based on observations. 

For example, the Langevin type models are stochastic differential equations de- 
scribing various phenomena in physics, biology, and other flelds. SDEs are used to 
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model various price processes, exchange rates, and interest rates, among others, in 
finance. Noises in these SDEs may be modeled as a generalized time derivative of 
some distinguished stochastic processes, such as Brownian motion (BM) or other 
processes. 

We are interested in estimating parameters contained in the stochastic differen- 
tial equation p , so that we obtain computational models useful for investigating 
complex dynamics under uncertainty. 

Theoretical results on parameter estimations for SDEs driven by Brownian mo- 
tion are relatively well developed,Sl2lIll2lIl2l33l46] ^^^^ 

various numerical simulations 
for these parameter estimationJ^^^EiEsl are also implemented. Sed^ for a more re- 
cent review about estimating and computing uncertain parameters, when dynamical 
systems are submit to colored or non-Gaussian noises. 

These research works on estimating uncertain parameters in dynamical systems 
are based on observations of sample paths. In the next section, we devise a method 
to estimate uncertain parameters based on observations of probability distributions 
of the system evolution. 



3 . Quantifying uncertainty in the space of probability measures 

Consider a dynamical system with model uncertainty, modeled by a scalar SDE 

dX = b{X)dt + a{X)dBt, X{0)=xo, (3.1) 

where the drift b{X) and diffusion a{X) contain uncertain parameters, to be esti- 
mated based on observations of probability distributions (i.e., probability measures) 
of the system paths Xt- 

To this end, we need to introduce the Hellingcr distanc^^^l between two prob- 
ability measures. It is used to quantify the similarity between two probability 
distributions. 

This is a metric in the space of probability measures. 



3 .1. Stationary case 

Under appropriate conditions on b and a (see,'^p.l70.), such as, 6 < and a ^ as 
well as some smoothness requirements, there exists a stationary probability density 
p{x) for the SDE (|3 .ip . as a solution of the steady Fokker-Planck equation, 

where the positive normalization constant C is chosen so that p > and J^p{x)dx = 
1, i.e., 

^^V/ , dx. 
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Note that xq here may be an arbitrary point so that the integral /j^^ ^2(^j dy exists 
(say, take = if that is convenient). 

Example 3 .1. (a) A special case: Langevin equation 

dX = ~bXdt + dBt, 

with parameter b > 0. Given an "observation" of the stationary probability density 
q[x) = -^e~^ • Find a 6 so that the Hellinger distance F{b) = \ /r[\/p(^ — 
^J q(x)Ydx is minimized. 

Hint: Insert p^q into the Hellinger distance. F{b) is an algebraic function of 
parameter & > and thus we use deterministic calculus to find a minimizer b 
(possibly by hand, or Matlab if needed). Note: /jje"^ dz = ^/tt. 

(b) A more general case: 

dX = b{X)dt + dBt, 

with function b{x) < 0. Given an "observation" of the stationary probability density 
q{x) = ^ ^^^2 (the Cauchy distribution). Find a function b{x) < so that the 
Hellinger distance F{b{x)) = \ j^\J p{x) — \J q{x)\'^dx is minimized. 

Hint: Insert p, q into the Hellinger distance. F{b{x)) is a functional of b{x) and 
thus we use calculus of variations (on F{b{x))) to find a minimizer b. Derive the 
Euler-Lagrange equation to be satisfied by b{x), together with appropriate boundary 
conditions (needed for p > and J^pdx = 1). Can you devise an algorithm to 
simulate for b{x)7 

Solution: 

(a) From the Fokker-Planck equation of 

dXt = -bXtdt + dBt, X() = 0, 
we have the stationary probability density 

To minimize the Hellinger distance F{b), we need 

F'(6) = - / e ( — = —x'^)dx -= e ^ {-:b ^ - —b-i)dx 

= i=6i(l + 6)-i-^rl(l + 5)-i=0. 
Therefore, 
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(b) From the Fokkcr-Planck equation of 

dXt = b{Xt)dt + dBt, Xq = 0, 
we could have the stationary probabihty density 

J —OD 

Using calculus of variation, we need 



F{b{x) + eh{x)) 



^2j^{b[v)+eh{y))dy 2efo Wv^+^f^iv^^dx 



0. 

This is equivalent to 

S e^S "Md^ _ V-^^ — . / (e^fo / )dz = 0. 

\ \/7r(l + a;2) I Jr\ Jz ) 

When changing the order of integration (first on y and z then on y and x), we 
have 



dec 



holds for all h{y). 
Therefore, 



J/ 



And since H e^/o K^)d^dz > 0, we have 



oo 



b{w)dw I g/; b{w)dw _ V-i«- \dx = 

for y G R. When taking the derivative of the above equation with respect to y, we 
have 
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Therefore, after taking 'square' and ' In' on both sides of the above equation, we 
have 

Jo Jr 
Taking derivative with respect to y, we have 

7r(l + 

only need b{w)dw < for all y e R. 



3 .2. Time- dependent case 
Consider a scalar SDE 

dX = b{X)dt + a{X)dBt, X(0) = X. (3 .2) 

The Fokker-Planck equation 

Pt = ]^{'^'^{x)p{x,t))^^ - {b{x)p{x,t))x, (3 .3) 

Example 3 .2. Consider a scalar SDE 

dX = -b sm{X)dt + V2dBt, X{0) = 0. (3 .4) 

(1) Assumc that an observation obtained for p to be 

gi(x,<)=.-l=e-^. (3.5) 

(2) Assunic that another observation obtained for p to be 

n{t + x'^)' 

The Fokker-Planck equation for p .4p is 

Pt = Pxx + {b ain{x)p)^. (3 .7) 

Find b based on the above two observation. 



Solution: In this case, we define the Hellinger distance: 

/oo 
{y^qt{x,t) - \/p{x,t,b)f'dx 
-oo 
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where i = 1,2. We need to calculate b by minimizing Hi(h). 

Now it's not easy to get the analytical solution oip{x, t), so we do some computer 
simulations instead. Assuming that the initial distribution p(a;, 0) is a delta function 
do on (—00,00), we can solve the above Fokker-Planck equation numerically, the 
solution is in the following figures. 





Fig. 0.1. Solutions of Fokkor-Planck equation when b 
30(right). 



1. In the plot, t = 5 (left) and t 



For the observation qi{x,t), it's the solution of pt = Pxx with delta function 
5q as the initial condition. Hence, the corresponding SDE of qi{x,t) is dX = 
dBt, X{0) = 0. We have its distribution shown below: 
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Fig. 0.2. Observation distribution qi(x,t). In the plot, t = 5 (left) and t = 30(right). 



By the definition Hi{b) — maxtgjo x] /,^(-\/ 9i(a;, — ■\/p(a;, t, b))'^dx^ we have 
the plot of Hi{h) in the following figures. And whatever T is , Hi{b) is always 
minimized when b = 0. 

For the observation q2{x,t), we have its distribution like 

By the definition H2{b) — maxjgjQ j-] J^oo^\/Q2{x, t) — \/p{x^ t, b))'^dx, we have 
the plot of H2{b) in the following figures. So we can see that if T = 5, H2{b) is 
minimized when & = 0.7 and if T = 30, H2{b) is minimized when b = 0.6. 
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